Universal Scaling of Wave Propagation Failure in Arrays of Coupled Nonlinear Cells 
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We study the onset of the propagation failure of wave fronts in systems of coupled cells. We 
introduce a new method to analyze the scaling of the critical external field at which fronts cease 
to propagate, as a function of intercellular coupling. We find the universal scaling of the field 
throughout the range of couplings, and show that the field becomes exponentially small for large 
couplings. Our method is generic and applicable to a wide class of cellular dynamics in chemical, 
biological, and engineering systems. We confirm our results by direct numerical simulations. 



PACS: 87.18.Pj, 82.40.Bj, 87.19.Hh, 84.30.-r 

The impact of discreteness on the propagation of phase 
fronts in biophysical, chemical, and engineering systems 
has been intensively studied during the last decade. 
Among the diverse examples are calcium release waves in 
living cells m, reaction fronts in chains of coupled chemi- 
cal reactors ||2]J^ , arrays of coupled diode resonators , 
and discontinuous propagation of action potential in car- 
diac tissue All these disparate systems share a 
common phenomenon of wave front propagation failure, 
independently of specific details of each system. Re- 
cently this effect has drawn considerable attention (see, 
e.g., Refs. Numerous experimental evidences show 
that the propagation failure occurs at finite values of the 
coupling strength (a critical coupling) [^^|J^ • This is con- 
trary to continuous systems, where wave fronts propagate 
for arbitrary couplings pc| ]. A challenging problem is to 
establish the universal properties of the critical coupling; 
this is crucial for making predictions of qualitatively dif- 
ferent regimes of system dynamics. 

In this Letter we consider the universal behavior of 
phase separation fronts in one-dimensional nonlinear dis- 
crete systems in an external field. We study the propa- 
gation failure transition for a class of simple dynamical 
models describing experimental observations in arrays of 
coupled nonlinear cells, such as chains of bistable chemi- 
cal reactors , systems of cardiac cells , etc. A new 
analytical method is presented to study generic proper- 
ties of the critical external field of the transition. We 
find, using this method, how the critical field scales with 
the intrachain coupling. This method is applicable for 
a wide range of the couplings. We confirm our analyti- 
cal predictions by direct numerical simulations of the full 
system. Our model in general is given by the following 
set of coupled nonlinear equations: 
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Here w„ is the order parameter at the n-th site, 7 is 
the damping coefficient, C is the coupling constant, and 
G{un,E) is the onsite potential, where E is the applied 
field. The potential has at least two minima separated by 



a barrier ub ■ The external field E is responsible for the 
energy difference between the minima. This provides for 
one globally stable and one metastablc minimum, u = u+ 
and u = , respectively. Phase fronts connect the two 
minima and tend to propagate, to increase the size of the 
energetically more favorable phase, w+ . 




FIG. 1. Stationary front near the propagation failure tran- 
sition in Eqs. (^). Parameters are (3 = 0.19 , E = 0.4 . Solid 
line represents the front solution. Dashed lines correspond to 
the minima ii+ and and the barrier ub ■ 

The mechanism of front propagation is manifested in 
the competition between the system discreteness and the 
driving field E . This competition gives rise to the propa- 
gation failure at the critical field Ec , which depends upon 
the intercellular coupling . For E > Ec , the front prop- 
agates at the velocity V (see, e.g., Ref. ||ll|) vanishing at 
the transition point, E ^ E^ . 

The results presented in this Letter are generic and 
apply to a wide class of systems, independently of the 
details, as long as the potential G{u,E) has the bistable 
structure described above [|l2|. We start our analysis 
with the sine-Gordon potential, G — K{1 — cos u — Eu) , 
with the factor K being the potential amplitude. This 
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potential was chosen in order to address the systems with 
phase dynamics, where the order parameter possesses a 
natural periodicity, such as arrays of Josephson junc- 
tions [ p^ . The dimensionless dynamics in this case are 
given by 



dun , 

= [i[Un+l + Un-l 



2Uri 



E — sin Un , (2) 



where time is rescaled as t — > tK/j and the dimensionless 
couphng is /3 = C/K . Our consideration is focused on 
the elementary fronts connecting two nearest minima out 
of the infinite set of the potential minima. The system 
is invariant with respect to the shift u — > w + 27r , so we 
choose without loss of generality u- = arcsini?, u+ = 
U- + 2n , ub — — arcsini? . This makes the dynamics 
effectively bistable. 

We introduce a novel analysis to find the universal de- 
pendence of the critical field Ec on the coupling (3 . For 
not too large /3 , the results are in a good quantitative 
agreement with the corresponding asymptotic descrip- 
tion. Finally we address a bistable fourth degree poly- 
nomial potential with the corresponding force in Eq. (|^) , 
-dG/du = -u{u - l)(u - 1/2 + E) . Such a poten- 
tial is applied, e.g., to describe the propagation failure 
in arrays of chemical reactors |^.|| and coupled diode 
resonators 1^,1). This illustrates the apphcability of our 
approach to generic potentials. The obtained analytical 
results for both potentials are justified by the numerical 
simulations of the full systems. 

For large dimensionless coupling, /3 1 , Eqs. (||) 
approach the continuous regime described by the over- 
damped sine-Gordon equation, ut = I3uxx — sin u + E . 
Here ec is a spatial coordinate standing for the contin- 
uous site number. This equation possesses a front so- 
lution, which propagates at nonzero velocity, for any 
driving field E . For small field, the front has the form 
u ~ 4arctan [exp(z/\/^)] , with the traveling wave coor- 
dinate z = x-Vt and velocity V ~ E^/P 0. Therefore, 
the critical field E^ vanishes in the continuous regime. 
For finite (3 , however, the external driving can be bal- 
anced by the effect of discreteness, thus leading to a prop- 
agation failure at a finite value oi E = Ec{/3) . Note that 
for a wide class of bistable systems, there is a global in- 
stability at E = Egi , above which the potential minima 
u± (or one of them) cease to exist. Typically the prop- 
agation failure occurs for Ec < Egi . In particular, for 
Eqs. @) one has Egi — 1 . 

In order to study the onset of the propagation failure, 
one has to consider the stationary case of Eqs. (^), i.e., 
dun/dt — (Fig |l|). We denote the site closest to the 
barrier separating two potential minima, as the "front 
site" and assign to it number n — Q . 

Before developing the main approach of the present 
study, we briefly present our results obtained by means 
of single-active-site theory, valid for the discrete regime 
of not too large coupling /3 . This theory is based on 



the fact that, for such (3, only the front site, n = 0, 
experiences nonlinearity, while all the other sites are close 
enough to either or m_ to be in a linear regime, see 
Fig 1^. Solving the linearized version of Eqs. (||) in the 
stationary case, one obtains u„ = A -I- Aexp(An) and 
M„ = 27r + A + _Bexp(— An) , for rt < and n > , 
respectively. Here A = arcsini?, A = arccosh[ 1/(2/3) -|- 
cos A ] . Matching these at n = leads to A = 27^ -\- B , 
which after substitution in the full nonlinear equation for 
Wo , yields 

- 2/? (7r + B)(e-^ - 1) +sin(B + A) =i; . (3) 

For E < Ec , Eq. (^ possesses two solutions, correspond- 
ing to stable and unstable fronts in the original problem 
(^. At the bifurcation point E = Ec these solutions 
merge and disappear, so that for E > Ec no station- 
ary solution to Eqs. (^ exists, and the front propagates. 
This implies that at the bifurcation point the derivative 
of the left hand side of Eq. (|^) with respect to B must 
equal zero. Then we obtain, after some calculations, an 
approximate expression for Ec : 



Ec 



V2V1 + 4/3-4/3-1 - /(/?) arccos [/(/3)] 



1 + fm 



(4) 



where /(/3) = 2 



1 - 



(3 + 1 



■ 2(3 + VTTTpi 

The graph of Ec{(3) given by (||) is shown in Fig. |^ 
(dotted line). We see from the figure that it agrees well 
with the results of numerical simulations of the full sys- 
tem (H), for small to moderate values of /3 (/3 ^ 0.5). This 
makes the single-active-site theory substantially more 
useful than a regular small (3 perturbation theory, which 
works only for much smaller (3 (^ 0.1). 

Although the single-active-site theory gives reasonable 
predictions for the propagation failure transition for mod- 
erate values of the intrachain coupling (3 , it does not 
provide a universal scaling of the critical field Ec with (3 
(see Fig. This motivates developing a general theory 
to describe the phenomenon of the propagation failure 
throughout the range of (3. Such a theory is presented 
below. 

A stationary front in Eq. (H) can be obtained as an 
extremum of the free energy 
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EUr, 



(5) 



For the case of zero field E, the front is always stationary, 
and no transition exists. In this case we have for the 
front, taking into account the effect of discreteness (see 
also H) 



w(n) — 



sin w(n) 
12/3 



(6) 



w{x) 



-2 arctan 



h sinh 



w{x) = 27r — w(— x) , 



X < , 
X > , 
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with the coefficient 6 = [1 - 1/(12/3) 

At zero field E , the free energy £ (||) possesses an infi- 
nite set of minima of equal depths, separated by barriers. 
Each of these minima corresponds to a stable stationary 
front and differs from the other fronts by a shift on an 
integer number of sites. Each barrier corresponds to an 
unstable front. When the external field E is applied, the 
set of minima is tilted, so that the difference between 
their depths is determined by the field value. As the 
field increases, each of the minima approaches its adja- 
cent barrier (down the energy landscape), until they fi- 
nally merge, at E — Ec ■ For E > Ec , the energy has no 
extrema left, and therefore the fronts propagate without 
limit. 

Our method is based on choosing profiles that 
can be parameterized by an effective coordinate a , and 
replacing the argument n of function w{n) in (H) with 
n + a . We show below that this choice allows us to 
find the critical field Ec . It can be demonstrated that, 
for E — 0, the points a = 1/2 + m (m is an integer) 
correspond to stable fronts, and a = to , to unstable 
fronts. The former and the latter are energy minima and 
maxima, respectively, satisfying d£ /da = . Note that 
the chosen parameterization by a may be interpreted as 
a continuous shift of a profile {un} ■ 

Our goal now is to evaluate the free energy landscape 
£(a) given by (|5|), for nonzero field E. When E = Ec 
the minima of £{a) merge with their adjacent maxima at 
the infiection points d'^E/da^ = 0. We have substituted 
the fronts (^ with w{n + a) into energy (||) and found Ec 
for various /9's, using Mathematica software. The result 
is shown in Fig. |^ (dashed line). 

In order to find the analytical form of £{a), we use the 
Poisson summation formula for infinite series: 



^F(n) 



Fix) 



1 



2 ^ cos(27rA:a;) 



fc=i 



dx , (7) 



where F(n) is the nth element of the series in (|). One 
finds that the integral J F{x)dx w Ea , up to a term 
independent of a. Then, after closing the integration 
contour in the complex plane, we see that fc > 1 terms in 
(0) are exponentially small, compared to the k — 1 term. 
This dominant term can be evaluated approximately to 
give the following expression for £{a) 



£{a) — 2'KEa + f7/?cos(27ra) exp 



(8) 



The factor f2 is a constant and can be found by compari- 
son of (H) with the results obtained by Mathematica (see 
Fig. ^) , which yields « 340 . Then, using the condi- 
tions d£ /da = d'^S / da^ = , we find 



1 

4 ' 



Ec = exp 



^li - 1/12 



(9) 



The general scaling Ec{P) of the propagation failure tran- 
sition given by (|^), virtually coincides with the one ob- 
tained by Mathematica, shown in Fig. ^ ]l6[ | . In the con- 
tinuous regime of large f3 the result (^) has the following 
simple asymptotic form Ec = f2/3exp[— tt^y^] . This im- 
plies that the critical field Ec decays exponentially for 
large couplings /3 . 

To confirm our general theory, we have performed nu- 
merical simulations of the full system (||). We have used 
the implicit second order integration method, in order to 
obtain the transition line Ec{P) ■ The results are pre- 
sented in Fig. 1^ (solid line). We see in the figure a good 
quantitative agreement of our analytical predictions with 
the results of the simulations. 
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FIG. 2. Critical curve of the propagation failure transi- 
tion for the sine-Gordon potential. Predictions of our gen- 
eral method and single-active-site theory are compared with 
numerical simulations of system (^). The system length is 
L = 100, and time step of the simulations dt = 0.0001. 

To address the propagation failure in arrays of chem- 
ical reactors [^]J^ and diode resonators [||J|], we turn to 
a fourth degree polynomial potential leading to a cubic 
force in Eqs. (0): -dG/du = -u{u - l){u - 1/2 + E) . 
The globally stable and metastable minima for this po- 
tential are u+ — 1 and m_ = , respectively; the bar- 
rier ub = 1/2 — E . Applying our general theory de- 
veloped above, we find the propagation failure transition 
line Ec{f3) . In Fig. |3| we show the results obtained by 
Mathematica (dashed line). We see again that the theory 
is confirmed by numerical simulations of the full system 
(solid line). The analytic expression of EdP) analogous 
to (^) is cumbersome so in this Letter we only give 
its asymptotic form for large (3 : Ec = il(3 exp{—T]^/J]) . 
This form has the same structure as the asymptotics of 
critical curve (^ for sine-Gordon potential, which pro- 
vides a persuasive argument that the proposed approach 
is generic, independent of the specific potential. The val- 
ues of H. and 77 are found to be « 429 , 77 w 21.9 . 
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FIG. 3. Critical curve of the propagation failure transition 
for the fourth degree polynomial potential. The system length 
and and time step of the simulations are the same as in Fig. bl 



To conclude, we have introduced a general method to 
analyze the transition of front propagation failure in ar- 
rays of coupled nonlinear cells. Using this method we 
have determined, for bistable dynamics of the cells, how 
the critical external field of the transition scales with in- 
trachain coupling. To demonstrate the generality of the 
new method, we have carried out our analysis for two dif- 
ferent effectively bistable potentials: (i) the sine-Gordon 
potential, useful e.g., for describing arrays of Josephson 
junctions [p^ ; and (ii) a fourth degree polynomial po- 
tential applicable to chains of chemical reactors and 
coupled diode resonators Our theoretical predic- 

tions have been confirmed by numerical simulations of 
the full systems. After this study of the propagation 
failure problem in the one-component systems (|l|), the 
following question arises: what is the mechanism of this 
phenomenon in systems described by more complicated 
dynamics? For example, the "single-pool-model" for in- 
tracellular calcium waves [^8| is represented by two dy- 
namic components: the calcium concentration and the 
fraction of vacant ionic channels. An important and vir- 
tually unexplored issue is how the effect of discreteness 
of ionic channels determines the onset of the propaga- 
tion failure of calcium waves. We expect that our tech- 
niques can provide new insights into this and other types 
of more complicated chains of coupled nonlinear elements 
in biophysics, chemistry, and engineering. A further chal- 
lenging direction of study here is to extend our analysis 
of the propagation failure to two- and three-dimensional 
systems. This would allow one to understand the effect of 
discreteness on the dynamics and stability of such struc- 
tures as spiral and scroll waves in cardiac tissue |^,|l9|,0 . 
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